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ABSTRACT 

Much of the dynamical structure of the Kuiper belt can be explained if Nep¬ 
tune migrated over several AU, and/or if Neptune was scattered to an eccentric 
orbit during planetary instability. An outstanding problem with the existing for¬ 
mation models is that the distribution of orbital inclinations they predicted is 
narrower than the one inferred from observations. Here we perform numerical 
simulations of Kuiper belt formation starting from an initial state with Neptune 
at 20 < aN,o < 30 AU and a dynamically cold outer disk extending from be¬ 
yond aN,o to 30 AU. Neptune’s orbit is migrated into the disk on an e-folding 
timescale 1 < r < 100 Myr. A small fraction (~10“^) of the disk planetesimals 
become implanted into the Kuiper belt in the simulations. By analyzing the 
orbital distribution of the implanted bodies in different cases we hnd that the in¬ 
clination constraint implies that r > 10 Myr and aN,o ^ 25 AU. The models with 
r < 10 Myr do not satisfy the inclination constraint, because there is not enough 
time for various dynamical processes to raise inclinations. The slow migration 
of Neptune is consistent with other Kuiper belt constraints, and with recently 
developed models of planetary instability/migration. Neptune’s eccentricity and 
inclination are never large in these models (cn < 0.1, zn < 2°), as required to 
avoid excessive orbital excitation in the >40 AU region, where the Cold Classicals 
presumably formed. 


1. Background 

The Kuiper belt is a diverse population of trans-Neptunian bodies (Figured]). Based on 
dynamical considerations, the Kuiper Belt Objects (KBOs) are classihed into several groups: 
the resonant populations, classical belt, scattered/scattering disk, and detached objects (also 
known as the fossilized scattered disk). See Gladman et ah (2008) for a formal dehnition 
of these groups. The resonant populations are a fascinating feature of the Kuiper belt. 


- 2 - 


They give the Kuiper belt an appearance of a bar code with individual bars centered at 
the resonant orbital periods. While Pluto and Plutinos in the 3:2 resonance with Neptune 
(orbital period ~250 years) are the largest and the best-characterized resonant group, nearly 
every resonance hosts a large population of bodies. The resonant bodies are long-lived, 
because they are phase-protected by the resonance from close encounters with Neptune. 
The orbits of the scattered/scattering disk objects, on the other hand, evolved and keep 
evolving by close encounters with Neptune. These objects tend to have long orbital periods 
and to be detected near their orbital perihelion when the heliocentric distance is ~30 AU. 
Their neighbors, the detached objects, have a slightly larger perihelion distance than the 
scattered/scattering objects and semimajor axes beyond the 2:1 resonance (a > 47.8 AU). 
The detached objects probably suffered close encounters with Neptune in the past, were 
scattered to orbits with large semimajor axes and eccentricities (e > 0.24 dehnes them in 
Gladman et al. 2008), but then they became “detached” from Neptune when some process 
increased their perihelion distance (or when Neptune’s orbit circularized; Levison et al. 
2008). 

The classical belt is a population of trans-Neptunian bodies dynamically defined as 
having non-resonant orbits with perihelion distances that are large enough to avoid close 
encounters with Neptune. They can be thought as of being related to the detached objects 
but having orbits with modest orbital eccentricities (e < 0.24 according to Gladman et al. 
2008). Here we consider the main classical belt located between the 3:2 and 2:1 resonances 
with Neptune (39.4 < a < 47.8 AU), because this is where most known classical objects 
reside. It is useful to divide the main belt into the dynamically “cold” and “hot” components, 
mainly because the inclination distribution of the main belt orbits is bimodal (Brown 2001), 
hinting at different dynamical origins for these components. Here we adopt an approximate 
separation, with Gold Glassicals (GGs) being defined as having i < 5° and Hot Glassicals 
(HGs) as i > 5°. Note that this dehnition is somewhat arbitrary, because the continuous 
inclination distribution near i = 5° indicates that mixing between the two components must 
have happened (e.g., Morbidelli et al. 2008, Volk & Malhotra 2011). 

While the HGs share many similarities with other dynamical classes of KBOs (e.g., 
scattered disk, Plutinos), the GGs have several unique properties. Specihcally, (1) the GGs 
have distinctly red colors (e.g., Tegler & Romanishin 2000) that may have resulted from 
space weathering of surface ices, such as ammonia (Brown et al. 2011), that are stable 
beyond ~35 AU. (2) A large fraction of the 100-km-class GGs are wide binaries with nearly 
equal size components (Noll et al. 2008a,b). (3) The albedos of the GGs are generally higher 
than those of the HGs (Brucker et al. 2009). And finally, (4) the size distribution of the 
GGs is markedly different from those of the hot and scattered populations, in that it shows 
a very steep slope at large sizes (e.g., Bernstein et al. 2004, Fraser et al. 2014), and lacks 
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very large objects (Levison & Stern 2001). The most straightforward interpretation of these 
properties is that the CCs formed and/or dynamically evolved by different processes than 
other trans-Neptunian populations. 

Following the pioneering work of Malhotra (1993, 1995), studies of Kuiper belt dynamics 
first considered the effects of outward migration of Neptune that can explain the prominent 
populations of KBOs in major resonances (Hahn & Malhotra 1999, 2005; Chiang & Jordan 
2002; Chiang et al. 2003; Levison & Morbidelli 2003; Gomes 2003; Murray-Clay & Chiang 
2005, 2006). With the advent of the notion that the early solar system may have suffered 
a dynamical instability (Thommes et al. 1999, Tsiganis et al. 2005), the focus broadened, 
with the more recent theories invoking a transient phase with an eccentric orbit of Neptune 
(Levison et al. 2008, Morbidelli et al. 2008, Batygin et al. 2011, Wolff et al. 2012, Dawson 
& Murray-Clay 2012). 

The emerging consensus is that the HCs, together with the resonant, scattered and 
detached populations, formed in a massive planetesimal disk at <30 AU, and were dynami¬ 
cally scattered onto their current orbits by migrating (and possibly eccentric) Neptune, while 
the CCs formed at >40 AU and survived Neptune’s early “wild days” relatively unharmed 
(Batygin et al. 2011, Wolff et al. 2012). The main support for this model comes from 
the unique properties of the CCs, which would be difficult to explain if the HCs and CCs 
had similar formation locations (and dynamical histories). For example, the wide binaries 
observed among the CCs would not survive scattering encounters with Neptune (Parker & 
Kavelaars 2010). Moreover, if the CCs evolved from the high-eccentricity Neptune-crossing 
orbits, this process should produce a gradient in e with more orbits having large e and fewer 
orbits having small e. The CCs do not show such a trend. Instead, low eccentricities prevail 
in that population. 


2. The Inclination Problem 

The inclination distribution of various populations in the Kuiper belt can be represented 
by N{i) dt = siniexp(—i^/2(T/) dt, where cij is a parameter (Brown 2001). In the main belt, 
the inclination distribution is bimodal and two components are needed: ai ~ 2.0° for the low- 
i CCs and = 8°-17° for the high-i HCs (Brown 2001, Kavelaars et al. 2008, 2009, Gulbis 
et al. 2010). The low inclinations of CCs are in line with the expectation that they formed 
from a dynamically cold disk at >40 AU, and their orbits were never excited too much by 
subsequent dynamical processes. The high inclinations of the HCs, on the other hand, are 
more challenging to explain (see below). Moreover, there is some evidence from high-latitude 
surveys that siniexp(—i^/2(T/) may be somewhat inadequate, because the drop-off at large 
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values of i is probably steeper than expected from this functional dependence (Petit et al. 
2015). For this reason, it is possible that N{i) di = sinzexp(—(i — ioY/erf) di with zq ^ 5° 
may better represent the underlying distribution. 

The HC distribution with relatively high orbital inclinations is shared among several 
other Kuiper belt populations as well, including Neptune Trojans (NTs), and the resonant 
and scattered objects. Eight NTs are currently known. Four of them have orbital inclinations 
i < 10°, and four have 25 < i < 30°. This could mean that the distribution is bimodal, but 
Parker (2015) showed that the bimodality of the underlying inclination distribution cannot 
be demonstrated with conhdence from the existing data. If the distribution is parametrized 
by a single term, N{i)di = smiexp{—i‘^/2af)di, the NTs are inferred to have exj > 11° 
with a 95% conhdence (Parker 2015). Plutinos in the 3:2 resonance with Neptune are well 
represented by a single term with exj ~ 11° according to Gulbis et al. (2010), or exj ~ 15° 
according to Kavelaars et al. (2008) and Gladman et al. (2012). Interestingly, the GG-like 
component with low orbital inclinations is not found in the Plutino population. 

The wide inclination distribution of the HGs, NTs, Plutinos and other resonant pop¬ 
ulations poses an important constraint on dynamical models of Kuiper belt formation ^ It 
implies that some dynamical process must have increased the inclinations by 10-15° on aver¬ 
age, and by ~30° at least in some cases. For example, the inclination constraint can be used 
to rule out a model in which these populations arise from a dynamically cold planetesimal 
disk at >30 AU, simply because the orbital inclinations are not excited in this region dur¬ 
ing Neptune’s migration (passing mean motion resonances do not affect inclinations much). 
Hahn & Malhotra (2005) investigated this issue in detail and found that starting with mod¬ 
erately excited orbits of planetesimals (e.g., i < 10°) does not resolve the problem, because 
the hnal inclination distribution is still not wide enough. The only way to make things work 
in the context of the Hahn & Malhotra model would be to assume that the inclinations were 
already high before Neptune’s migration, but that seems unsatisfactory, because it is not 
clear how the inclinations could have been excited beforehand. Lykawka & Mukai (2008) 
considered dynamical effects of an additional planet in the trans-Neptunian region. They 
found that this putative planet could help to excite inclinations, but the inclination distri¬ 
bution obtained in their model was not wide enough to match observations well. It lacked 
orbits in the classical belt with i > 15°, while these orbits are in fact common. 

Given the difficulties described above, various theoretical models considered the forma¬ 
tion of the HGs and resonant populations from a massive planetesimal disk at <30 AU (the 


^Assuming the current configuration of planets, long-term orbital dynamics in the Kuiper belt region 
cannot explain the high inclinations of the KBOs (Kuchner et al. 2002, Li et al. 2014). 
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outer edge of the massive disk is constrained to ~30 AU by Neptune’s present orbit; Gomes 
et al. 2004). To reach >40 AU, planetesimals must be radially displaced. Levison & Mor- 
bidelli (2003) considered a scenario in which objects were pushed out by the 2:1 resonance 
with Neptune. This could work only if Neptune was initially inside ~19 AU, such that the 
2:1 resonance fell inside the massive disk’s outer boundary at ~30 AU. Gomes (2003), on the 
other hand, suggested that bodies were hrst scattered to >40 AU by having close encounters 
with Neptune, and became dynamically decoupled from Neptune while Neptune was still 
migrating. If so, the HGs and resonant populations would be close relatives of the scattered 
disk objects. The exact nature of the decoupling process for the HGs is uncertain (Gomes 
2003, Levison et al. 2008, Dawson & Murray-Glay 2012), but recent work suggests that 
capture into mean motion resonances with Neptune (2:1, 5:3, 7:4, etc.), and secular/Kozai 
cycles inside the resonances may have played an important role (Brasil et al. 2014b). 

These results could help to resolve the inclination problem discussed above, because the 
orbital inclinations can be excited when bodies undergo a series of close encounters with 
Neptune (Gomes 2003). In addition, inclinations are increased for orbits that suffer Kozai 
cycles, because to decouple from Neptune, the eccentricity must drop, and the inclination 
would therefore rise (due to the anticorrelated behavior of e and i caused by the Kozai cycles; 
Kozai 1962). It remains to be shown, however, how these processes operated to affect the 
Kuiper belt, and how the early evolution of planetary orbits is constrained by the dynamical 
structure of the Kuiper belt. 

As demonstrated in Figure [H the inclination problem is not resolved by simply postu¬ 
lating that much of the Kuiper belt has been implanted from <30 AU (e.g., see discussion 
in Petit et al. 2011). Here we choose to illustrate the inclination problem with Plutinos, 
because the 3:2 resonance population is characterized much better than any other resonant 
population (e.g., Gladman et al. 2012). Also, Plutinos do not show the bimodal inclination 
distribution of the classical belt, so we do not need to worry about the overlap of different 
groups. The parameters of the numerical model from Figure [2] are similar to those used 
in Levison et al. (2008). We used a fast migration regime with an exponential e-folding 
timescale r = 1 Myr (see next section for our model description) to illustrate that such a 
fast migration of Neptune leads to an implausible result. There is simply not enough time 
in this case to substantially raise the orbital inclinations. 

Here we perform numerical simulations to investigate the inclination problem in detail. 
Our method and constraints are described in Sections 3 and 4, and the results are presented 
in Section 5. We hnd that bodies starting with a < 30 AU can be implanted into the 
Kuiper belt by hrst being scattered by Neptune to > 30 AU, and then decoupling from 
Neptune by various resonant effects while Neptune is still migrating (see Section 5.2). We 
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show that the inclination constraint implies a prolonged phase dnring which Neptnne slowly 
migrated (e-folding migration timescale t > 10 Myr) before reaching its cnrrent orbit at 
a = 30.1 AU@ The main effect of slow migration is that orbits are allowed to deconple 
from Neptnne relatively late dnring the migration process. Conseqnently, Neptnne is given 
more time to act, via scattering enconnters, on the sonrce popnlation, thns increasing the 
orbital inclinations of bodies before they are implanted into the Kniper belt. The model with 
Neptnne’s slow migration is also consistent with other Kniper belt constraints (Section 5). 
Varions implications of this resnlt are discnssed in Section 6. 


3. The Integration Method 

Onr nnmerical integrations track the orbits of fonr planets (Jnpiter to Neptnne) and 
a large nnmber of test particles representing the onter planetesimal disk. To set np an 
integration, Jnpiter, Satnrn and Uranns were placed on their cnrrent orbits]^ Neptnne was 
placed on an orbit with semimajor axis aN,o; eccentricity eN,o; and inclination q- To cover 
the parameter space, we set aN,o = 22, 24, 26 or 28 AU, eN,o = 0, 0.1 or 0.3, and z = 0° or 
5°. We tested many different combinations of these parameters to nnderstand their role in 
Kniper belt formation. The cases with aN,o = 22 AU or 24 AU, eN,o = 0 and i = 0° are a 
good proxy for the initial conditions of Hahn & Malhotra (2005), who stndied a long-range 
migration of Neptnne, and for the instability/migration models developed in Nesvorny & 
Morbidelli (2012). The case with aN,o = 28 AU and eN,o = 0.3 is similar to rnns A, B and 
C in Levison et ah (2008), who motivated their choice by the strong planetary instability 
occnrring in the Nice model (Tsiganis et ah 2005, Morbidelli et ah 2007). The cases with 
eN,o ^0.1 were favored by Dawson & Mnrray-Clay (2012) from the CC-related constraints. 

The swift_rmvs4 code (Levison & Dnncan 1994) was nsed to follow the evolntion of 
planets and disk particles. The swift_rmvs4 code was modihed to inclnde hctitions forces 
that mimic Neptnne’s radial migration and damping. These forces were parametrized by 
exponential e-folding timescales Ta, Te and r*, where Tq controls the radial migration rate, 
and Te and r, control the damping rates of e and i. Here we set Ta ^ ^ Tt (= r), becanse 

snch ronghly comparable timescales were snggested by previons work. Specifically, we nsed 


correlation between the width of the inclination distribution and the migration timescale/time of 
capture were previously reported by Malhotra (1995) and Levison et al. (2009). 

•^The dependence of the results on the orbital behavior of Jupiter, Saturn and Uranus was found to be 
minor. We determined this by comparing our nominal results with fixed orbits to those obtained when these 
planets were forced to radially migrate. 
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r = 1, 3, 10, 30 and 100 Myr, where r = 1 Myr corresponds to the case considered by Levison 
et ah (2008), while t > 10 Myr is preferred from the instability simnlations of Nesvorny 
Sz Morbidelli (2012). By hne tnning the migration parameters, the hnal semimajor axis of 
Neptnne was set to be within 0.05 AU of its cnrrent mean on = 30.11 AU, and the orbital 
period ratio, Pn/-Pu, where Pn and Pu are the orbital periods of Neptnne and Uranus, was 
adjusted to end up within 0.5% of its current value (Pn/Pu = 1-96). 

Each simulation included one million disk particles distributed from just outside Nep¬ 
tune’s initial orbit to 30 AU. Their radial profile was set such that the disk surface density 
E oc 1/r, where r is the heliocentric distance. The large number of disk particles was needed 
because the capture probability in different parts of the Kuiper belt is expected to be ~10“^- 
10“^. With 10® disk particles initially, this yields ~100-1000 captured particles, and allows 
us to perform a detailed comparison of the model results with observations (Section 4). The 
disk particles were assumed to be massless such that their gravity does not interfere with 
the migration/damping routines. This means that the precession frequencies of planets are 
not affected by the disk in our simulations, while in reality they were. This is an important 
approximation, because the orbital precession of Neptune during its high-eccentricity phase 
can influence the degree of secular excitation of the CCs (Batygin et ah 2011). 

All simulations were run to 1 Gyr. The interesting cases were extended to 4 Gyr with 
the standard swift_rmvs4 code (i.e., without migration/damping in the 1 to 4 Gyr interval). 
We performed eighteen simulations in total. Three of these runs were designed to test the 
reproducibility of the results. While the results concerning the orbital distribution of bodies 
implanted into the Kuiper belt (Sections 5.1-5.3) were found to be strictly reproducible, 
the efficiency of capture in the Kuiper belt can vary by a factor of a few depending on the 
behavior of Neptune’s eccentricity during the simulation (see discussion in Section 5.4). 

An additional uncertain parameter concerns the dynamical structure of the original 
planetesimal disk. It is typically assumed that the disk was dynamically cold with orbital 
eccentricities e < 0.1 and orbital inclinations i < 10°. Some dynamical excitation could have 
been supplied by scattering off of Pluto-sized and larger objects that presumably formed 
in the disk (Stern & Golwell 1997, Kenyon et ah 2008)0 The magnitude of the initial 
excitation is uncertain, because it depends on several unknown parameters (e.g., the number 
of massive objects in the disk). Here we operate under the assumption that the orbital 
inclinations of disk particles were relatively small initially, and were excited during the main 


"‘The escape velocity from Pluto is 1.2 km s“‘, about 20% of the Keplerian orbital speed at 25 AU. 
Therefore, eccentricities up to ~0.2 and inclinations up to ~12° can be expected from a surface-grazing 
flyby near a Pluto-class object. 



stage of planetary instability/migration, when bodies were implanted into the Kniper belt. 
This is a reasonable assnmption, given that the notion of planetary instability/migration was 
developed, among other reasons, to explain the complex orbital strnctnre of the Kniper belt. 
It wonld thns seem nnsatisfactory to “resolve” the inclination problem discnssed in Section 
2 by postnlating that the inclinations were already large initially (nnless it is explained how 
that happened). The initial eccentricities and initial inclinations of disk particles in onr 
standard simnlations were distribnted according to the Rayleigh distribntion with = 0.05 
and Oi = 2°, where a is the usnal scale parameter of the Rayleigh distribntion (the mean 
of the Rayleigh distribntion is eqnal to For completeness, we also tested (Tg = 0.1 

and (Tg = 0.2, and (Tj = 5° and cxj = 10° in several cases. 


4. Constraints and the CFEPS Detection Simnlator 

The resnlts of the simnlations described in the previons section were compared to ob¬ 
servations. We paid special attention to the inclination problem described in Section 2, bnt 
also made snre that the best models identihed here are consistent with other Kniper belt 
constraints. Onr primary constraints were: 

1. The captnre efficiency and orbital distribntion of HCs. According to Fraser et ah 
(2014), the HCs contain a mass Mhc — 0.01 Me, where Me = 6 x 10^^ g is the 
Earth’s mass. With Mdisk = 20 Me, the captnre probability of HCs wonld therefore be 
Fhc — 0.01/20 = 5 X lO”"^. This estimate is probably at least a factor of ~2 nncertain, 
becanse both Mhc and M^isk are somewhat nncertain. The inclination distribntion of 
HCs obtained in the model is reqnired to be similar to the wide inclination distribntion 
inferred from observations. A detailed comparison is done with the Canada-France 
Ecliptic Plane Snrvey (CFEPS) detection simnlator (see below). The distribntion of 
a and e follows a trend seen in Fignre[31 where larger valnes of a correspond to larger 
valnes of e. This trend, which can be an important diagnostic of the implantation 
mechanism, mnst be reprodnced in a snccessfnl model. Also, the eccentricities of 
objects captnred in the main belt mnst reach below 0.05, as they do in reality. 

2. The CCs at 42 < a < 45 AU mnst snrvive and their orbits cannot be excited too 
mnch. Dawson & Mnrray-Clay (2012) snggested that the eccentricities of CCs were 
not excited above 0.05 in the inner part of the main belt (42 < a < 43.5 AU) and 
above 0.1 in the onter part (43.5 < a < 45 AU), becanse there appears to be a stable 
bnt nnpopnlated region above these limits. Morbidelli et al. (2014) demonstrated, 
however, that the 7:4 and 9:5 resonances conld have depleted the region in qnestion if 
Neptnne was on a somewhat eccentric orbit (cn — 0.1), when it reached = 28 AU, 
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and migrated slowly. Our results, discussed in Section 5, are in line with these hndings. 
We do not explicitly discuss the CCs in the following text, because our main results 
were obtained with cn < 0.1. We checked that the CCs are not excessively excited in 
this case, in agreement with Dawson & Murray-Clay (2012; also see Wolff et ah 2012). 

3. The capture efficiency and orbital distribution of the resonant objects. According to 
the CFEPS survey, there are ~3.5 as many HCs as Plutinos with absolute magnitude 
H < 8 (diameter D > 150 km for 0.05 albedo) (B. Gladman, personal communication; 
Petit et al. 2011, Gladman et al. 2012). This suggests a capture probability into the 3:2 
resonance of — 1-5 x 10“^ (estimate at least a factor of ~2 uncertain). The orbits 
of Plutinos show moderate to high inclinations, similar to those found for the HCs, 
and eccentricities mainly in the 0.1-0.35 interval. The populations in the 2:1 and 5:2 
resonances with Neptune are probably somewhat smaller (~2-4 times) than Plutinos 
(e.g., Gladman et al. 2014). The population of NTs is much smaller (Alexandersen et 
al. 2015), indicating a capture probability of the order of Pi-i ~ 10“®. All resonant 
populations have a wide inclination distribution (e.g., ai > 11° for NTs; Parker 2015). 

4. The existence and orbits of the detached objects. The detached objects have stable non¬ 
resonant orbits with semimajor axes beyond Neptune’s 2:1 mean motion resonance (a > 
47.8 AU) and perihelion distances up to g ~ 40 AU. These objects cannot be placed 
on their orbits in the current conhguration of the planetary orbits and thus provide 
an important constraint on any formation model. Levison et al. (2008) suggested 
that the detached disk was created during a phase when Neptune had a substantial 
orbital eccentricity (cn ~ 0.3) and was capable of scattering objects up to g ~ 40 AU. 
Here we show that the detached disk can be obtained even for cn — 0, assuming that 
r > 10 Myr (consistent with the condition required from the inclination constraint). 
The dynamical mechanism responsible for the formation of the detached disk (and 
HCs) is found to be a three-step process related to the capture of scattered bodies in 
migrating Neptune resonances (Section 5.2; Gomes 2003, Gomes et al. 2005, Brasil et 
al. 2014a). 

We used the CFEPS detection simulator (Kavelaars et al. 2009) to compare the orbital dis¬ 
tributions obtained in our simulations with observations. CFEPS is one of the largest Kuiper 
belt surveys with published characterization (currently 169 objects; Petit et al. 2011). The 
simulator was developed by the CFEPS team to aid the interpretation of their observations. 
Given intrinsic orbital and magnitude distributions, the CFEPS simulator returns a sample 
of objects that would have been detected by the survey, accounting for flux biases, pointing 
history, rate cuts and object leakage (Kavelaars et al. 2009). In the present work, we input 
our model populations in the simulator to compute the detection statistics. We then com- 
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pare the orbital distribution of the detected objects with the actual CFEPS detections using 
the Kolmogorov-Smirnov (K-S) test (Press et ah 1992). 

This is done as follows. The CFEPS simulator takes as an input: (1) the orbital ele¬ 
ment distribution from our numerical model, and (2) an assumed absolute magnitude {H) 
distribution. As for (1), the input orbital distribution was produced by a short integration 
starting from the hnal model state of the Kuiper belt. The orbital elements of each object 
were recorded at 100 yr intervals during this integration until the total number of recorded 
data points reached ~10^. Each data point was then treated as an independent observational 
target. We rotated the reference system such that the orbital phase of Neptune in each time 
frame corresponded to its ecliptic coordinates at the epoch of CFEPS observations. This 
procedure guaranteed that the sky positions of the objects in Neptune’s resonances were 
correctly distribnted relative to the pointing direction of the CFEPS frames. 

The magnitude distribution was taken from Fraser et ah (2014). It was assnmed to 
be described by a broken power law with N{H)dH = dH for H < Hb and 

N{H)dH = iQ‘^ 2 iH-Ho)+(ai-a 2 )(HB-Ho) ^ where ai and a 2 are the power-law 

slopes for objects brighter and fainter than the transition, or break magnitude Hb, and 
Hq is a normalization constant. Fraser et ah (2014) found that ai = 0.9, 0:2 = 0.2 and 
Hb = 8 for the HCs. In the context of a model where the HCs formed at < 30 AU, and 
were implanted into the Kniper belt by size-independent processes (onr integrations do not 
have any size-dependent component), the HC magnitnde distribntion should be shared by all 
populations that originated from < 30 AU (Morbidelli et al. 2009b, Fraser et ah 2014). We 
varied the parameters of the input magnitude distribution to nnderstand the sensitivity of 
the results to various assumptions. We found that small variations of ai, 0:2 and Hb within 
the nncertainties given in Fraser et al. (2014) have essentially no effect. 


5. Results 

5.1. A Reference Case 

We hrst discuss a reference simnlation with slow migration of Neptune (r = 30 Myr, 
on.o = 24 AU, eN,o = 0, iN,o = 0), to illustrate that the results of this model match the 
orbital structure of the Kuiper belt, including the wide inclination distribution of the HCs 
and resonant populations. Later, in Section 5.3, we will explain how the resnlts differ from 
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the reference case when varions model parameters, snch as r and aN ,05 variedjj FignreH] 
shows the orbital distribntion of the model orbits obtained with aN,o = 24 AU and r = 30 
Myr. This hgnre can be compared to Fignre [31 bnt note that some caution needs to be 
exercised in this comparison, because Figure [31 includes various observational biases, while 
Figure[3|does not. Also, the total number of points in the two plots is different (known KBOs 
with good orbits in Figure [3] and a fraction of the initial 10® disk particles in Figure [3|). 

The model results in Figure [H show a remarkable similarity to Figure [31 The orbital 
structure obtained in the model shows all main components of the Kuiper belt, including the 
resonant populations, classical belt, scattered and detached disks. The resonances such as 
the 5:4, 4:3, 5:3, 7:4 and 5:2 are also populated. The model orbits in the detached disk have 
perihelion distances reaching toward g ~ 40 AU, as they do in reality. [Note that the CCs 
are not shown in Figure [H because the model discussed here does not account for objects 
that formed beyond 30 AU (see discussion of the CCs in Section 1 and the description of the 
model in Section 3).] 

The model distribution of orbital inclinations in Figure [H covers the whole interval 
between 0 and 40°. This is a notable result because the original disk orbits had cxj = 2°. The 
orbital inclinations have therefore been signihcantly excited during the implantation process. 
The dynamical processes responsible for the implantation of objects in the Kuiper belt and 
their effects on the orbital inclination are discussed in Section 5.2. Here we hrst more carefully 
compare the model distribution with observations. To do this, the model distribution shown 
in Figure [H was passed through the CFEPS detection simulator. Figure [5] shows how the 
model detections compare with the actual CFEPS detections. The comparison is done 
separately for the 3:2 resonance and HCs. The reason for this is that the implantation process 
and stability properties can, and indeed do, produce differences between these populations. 
We do not show a similar comparison for other resonances, because the number of actual 
CFEPS detections in the 1:1, 2:1 and other resonances is very small (1 to 5 objects detected) 
and a rigorous comparison is therefore not possible at this time. 

Figure [5] shows that the orbital distribution of the detected objects in the 3:2 resonance 
agrees with the distribution of the CFEPS detections. The K-S test indicates that the 
model and observed distributions shown in Figures [5^ (eccentricity) and [Sb (inclination) 
have 67% and 84% probabilities, respectively, of being derived from the same underlying 
distribution. This is very good agreement]^ The model distribution of inclinations obtained 


^These results are strictly reproducible. There is enough information given in Section 3 for anyone to 
repeat our simulations and confirm the results. 

®To be more precise, a K-S probability of 0.67 (or 0.84) means the following. Assume that there was a 
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with r = 30 Myr is much wider than the one obtained for r = 1 Myr (Figure [2]), and matches 
observations very well. The eccentricity distributions are also very similar. Most Plutinos 
have e > 0.15 (>80% of detections). This characteristic is a consequence of the implantation 
mechanism, where orbits are deposited into the 3:2 resonance from the scattered disk and 
retain somewhat large eccentricities (see Section 5.2). Note that, with on.o = 24 AU, the 
3:2 resonance is initially outside the outer boundary of the planetesimal disk. Plutinos are 
therefore not captured in our model from the low-eccentricity orbits as in, for example, Hahn 
& Malhotra (2005). 

The agreement for the HCs is also good. The K-S test applied to the eccentricity 
and inclination distributions shown in Figures [51:; and[5li gives 37% and 75% probabilities, 
respectively. As for the inclination distribution of the HCs, our model with r = 30 Myr 
and aN,o = 24 AU predicts that about 10% of the CFEPS detections should have i > 30°, 
while no object was thus far detected by CFEPS with such a high inclination. This is not a 
problem, however, because the CFEPS detected only 10 HCs with i > 10°, and the statistical 
constraints for i > 30° are very weak. Por i < 10°, on the other hand, there are concerns with 
contamination from the CCs, which are not modeled here. The CCs should clearly have the 
dominant contribution to the statistics for i < 5°. For 5° < i < 10°, the situation is unclear. 
There are 11 CFEPS detections in this intermediate inclination range, which is similar to the 
number of detections for i > 10°. In our model with r = 30 Myr, however, these intermediate 
inclinations are not populated as much (we get about half of the expected detections). We 
believe this issue arises because our model with smooth migration of Neptune and a hxed 
value of r is only an approximation of the real evolution (see discussion in Section 6). 


5.2. The Implantation Mechanism 

We examined the orbital histories of test particles in the reference simulation and found 
that the implantation of bodies from < 30 AU into the Kuiper belt is in general a three-step 
process. The hrst two steps are common for the HCs and resonant populations; the third 
step is what distinguishes them. We first describe these steps and then illustrate them with 


single parent distribution for both the model and the observations. In particular, the model was a random 
sample containing, say, J entries, while the observations contained K entries. If we were to generate two 
random representations of this parent population, one with J entries and one with K, there would be a 67% 
(or 84%) chance that the comparison between these random populations would be worse than what we have 
in Figure Hh (or Figure [bh)- This holds despite the fact that these new distributions were directly derived 
from the parent. Therefore, the agreement between our model and the observations is very good. Indeed, 
any comparison with a K-S probability greater than ^0.1 should be considered acceptable. 
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a few examples. Specifically: 

Step 1: The disk planetesimals are scattered by Neptune from <30 AU to >30 AU. Their 
distribution resembles that of the scattered disk in that they populate the region with 
a > 30 AU and q < aN(t), where is the semimajor axis of migrating Neptune. 

Here we define the Intermediate Source Region or ISR as the orbital region with q < 
where QN(t) = aN(l + eN) is Neptune’s aphelion distance, and 40 < a < 47 AU. 
Most bodies implanted into the main belt evolved onto their present orbits via the ISR. 

Step 2 : The scattered bodies evolve onto orbits with large libration amplitudes in mean 
motion resonances. The secular dynamics inside the mean motion resonances is com¬ 
plex, including large-amplitude Kozai, apsidal and nodal cycles (e.g., Morbidelli et al. 
1995, Nesvorny & Roig 2000, 2001). These effects can act to decrease the orbital ec¬ 
centricity, thus decoupling the orbit from Neptune, on a characteristic timescale that 
is comparable to the period of the secular oscillations (>1 Myr). 

Step 3: If Neptune were not migrating, the evolution described in Step 2 would be reversible 
and bodies would be released, sooner or later, back to the scattered disk. With Nep¬ 
tune’s migration, however, two additional alternatives can happen: (1) the orbit can 
evolve to a smaller libration amplitude and stabilize inside the resonance, or (2) it 
can be released from the resonance with low eccentricity and can end up on a stable, 
HC-like orbit with q > 35 AU. 

The three-step implantation mechanism described above was originally proposed in Gomes 
Gomes et al. (2005). We will therefore call it the Gomes mechanism in the 

Figure [6] shows an example of a disk particle that was captured on a high-inclination 
and low-eccentricity orbit in the main belt. This case is a clear illustration of the three- 
step Gomes mechanism described above. Initially, the particle starts with a = 28.5 AU, 
e = 0.04 and i = 3.5°. It is scattered by Neptune and evolves into the scattered disk, where 
it remains until roughly t = 13 Myr after the start of the simulation. During this hrst stage, 
the eccentricity and inclination are excited by encounters with Neptune. The orbit is then 


(2003) and 
following]^ 


^The Gomes mechanism was previously shown to work in several self-consistent simulations of the plane¬ 
tary instability/migration, where the disk particles carried actual mass. Because a relatively small number of 
disk particles was used, however, Neptune’s migration was unrealistically grainy in these simulations. Here 
we show that the Gomes mechanism works, with a somewhat lower efficiency, even if Neptune’s migration 
was smooth. 
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captured in the 2:1 resonance and remains in the resonance until t = 22 Myr (see the libration 
of the resonant argument in panel f). Once captured, it undergoes Kozai oscillations (see 
the libration of the perihelion argument in panel e). The eccentricity drops from 0.4 to 0.06 
during this phase (13 to 22 Myr), while the inclination further increases from 18° to 28°. 
Finally, about 22 Myr after the start of the simulation, the particle is released from the 2:1 
resonance (see panels b and f) and lands on a main-belt orbit with a = 42.8 AU, e = 0.05 
and i = 28°. This orbit is stable for 4 Gyr. 

Figure [7] shows an example of a disk particle that was captured on a high-inclination 
orbit in the 3:2 resonance. In this case, the scattering phase lasted until about f = 65 
Myr. Both the eccentricity and inclination were strongly excited during this phase. At 
f ~ 65 Myr, the orbit entered into the 3:2 resonance. The resonant libration amplitude was 
initially variable but later evolved to become ~100° (panel f). The 3:2 resonant orbits with 
such amplitudes are stable (Nesvorny & Roig 2000). Indeed, the orbit stayed in the 3:2 
resonance for the whole duration of our integration (4 Gyr). During the capture into the 3:2 
resonance the orbit started showing Kozai cycles {u started librating in panel e while e and 
i in panels c and d exhibited signs of correlated oscillations). The Kozai oscillations with 
(full) amplitude of about 80° remained for the whole duration of the simulation. This case 
is reminiscent of Pluto, whose orbit also has Kozai cycles, but the hnal orbital inclination of 
this particle is considerably higher [i = 36° vs. Pluto’s i ~ 17°). 

Figure |H] is a case of Neptune Trojan. Akin to the example of the captured Plutino 
discussed above, the scattering phase lasts very long (~ 62 Myr). The eccentricity is excited 
but the inclination remains relatively low (<15°). Then, around f = 50 Myr, the orbit starts 
showing signs of the Kozai resonance (panel e) and undergoes two brief periods during which 
the 1:1 resonant angle librates (about t = 54 Myr and f = 60 Myr). The inclination rises and 
eccentricity drops in an anti-correlated pattern. The orbit is caught into the 1:1 resonance 
at t ~ 62 Myr, eventually stabilizes with a very small libration amplitude (~ 30°) around 
the leading Lagrangian point L 4 , and remains there for 4 Gyr. 

The three examples discussed above were selected from hundreds of similar cases to 
illustrate the Gomes mechanism for capture in the main belt and resonances, and the high 
inclinations that these orbits can reach, if the migration of Neptune is slow and the scattering 
phase lasts long. We found that the 2:1 resonance was typically involved in capture of orbits 
in the main belt, but resonances such as 7:4 or 9:5 also contributed. The orbits of the 
detached objects produced in our simulations follow a similar pattern, but the resonances 
responsible for raising their perihelia are different (e.g., 5:2, 7:3, 3:1). Moreover, relatively 
strong resonances such as the 3:1 are capable of producing the detached orbits with very 
large perihelion distances, and some of these orbits are reminiscent of that of 2004 XR 190 
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(Buffy) (a = 57.7 AU, q = 51.5 AU; Allen et al. 2006; see Figure HJ also Gomes 2011). 


5.3. The Inclination Distributions Obtained in Different Models 

The results of our simulations show that the inclination distribution of bodies implanted 
in the Kuiper belt depends both on r and aN,o; but the main effect is that of r. This is 
because the implantation of bodies in the Kuiper belt happens on a timescale comparable 
to r. With short r, the implantation must be fast and Neptune does not have much time to 
raise the inclinations of the scattered bodies. Consequently, the inclination distribution of the 
implanted bodies is narrow and clustered toward z ~ 0. If, on the other hand, r is long, the 
orbital inclinations of scattered objects can be substantially excited by Neptune encounters 
before these bodies are implanted into the Kuiper belt. The inclination distribution of the 
implanted bodies is thus wide in this case. Figures [2] and [5] illustrated the dependence on r 
for r = 1 Myr and r = 30 Myr. Figure IH] show this dependence for several additional cases. 

The inclination distribution of Plutinos obtained with different migration timescales is 
shown in Figure [9^. The results obtained for r < 10 Myr and any aN,o are clearly a poor £t, 
because they indicate cxj < 5°, while Plutinos have cxj > 10° (Kavelaars et al. 2008, Gulbis 
et al. 2010, Gladman et al. 2012). The best-fit distribution for the case with r = 10 Myr 
and aN,o = 24 AU is obtained for cxj ~ 10°. When this case is compared to the CFEPS via 
the CFEPS detection simulator iFigure fTObh we hnd that it can be ruled out with 99.6% 
conhdence. The results with r = 10 Myr and aN,o > 24 AU can be ruled out at even higher 
conhdence levels, because the inclination distributions of Plutinos obtained in those cases 
are narrower than the one obtained with r = 10 Myr and aN,o = 24 AU. We conclude that 
fast migration timescales with t < 10 Myr do not work for Plutinos. 

Longer timescales produce better results. The case with r = 30 Myr and on.o = 24 AU 
was discussed in Section 5.1 and clearly matches observations very well (Figure [5]). Interest¬ 
ingly, the inclination distribution obtained in the model also depends on aN,o- For example, 
the inclination distributions obtained with r = 30 Myr and aN,o = 26 AU or aN,o = 28 AU 
are much narrower (cTj ~ 10°) than the one obtained with r = 30 Myr and aN,o = 24 AU 
(cTj = 15 -20°). When compared to the CFEPS detections, these cases can be ruled out 
at a >99% conhdence level. This happens because bodies tend to be captured relatively 
early in these simulations when Neptune’s migration rate is still substantial. In contrast, the 
long-range migration with on.o ^ 25 AU offers more opportunity for capture at late times, 
and therefore leads to a wider inclination distribution that is more in line with observations. 
We therefore conclude from the inclination distribution of Plutinos that Neptune’s migration 
was slow (r > 30 Myr) and long range (aN,o ^ 25 AU). 
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An additional argument that favors a long migration timescale comes from the eccen¬ 
tricity distribution of Plutinos. With r < 10 Myr, the model eccentricity distribution is 
skewed toward large values of e when compared to observations. This can be demonstrated 
by comparing the model to the CFEPS detections via the CFEPS detection simulator. For 
example, for r = 10 Myr and aN,o = 24 AU iFigure fTOhh the K-S test applied to the ec¬ 
centricity distributions of the detected Plutinos gives only a 0.3% probability that the two 
distributions are the same. This mismatch is a consequence of step 2 of the capture process, 
discussed in Section 5.2, where there is not enough time available with rapid migration to 
decrease eccentricities. The eccentricities of captured objects therefore end up being too 
large. 

We now move to discussing the HC inclination distribution. The HCs are deposited 
into the main belt (40 < a < 47 AU) by mean motion resonances such as the 2:1 and 7:4. 
With aN,o = 28 AU, these resonances are located in the main belt with a = 44.4 AU and 
a = 40.7 AU, respectively. Therefore, in this case, the implantation of HCs into the main 
belt begins almost immediately after the start of a simulation (with a short delay required 
for Neptune to scatter bodies to the ISR; Section 5.2). The bodies implanted into the main 
belt during the initial stages will have small orbital inclinations and will skew the inclination 
distribution of HCs toward small values. Thus, even if long migration timescales are used in 
this case (e.g., r = 30 Myr or 100 Myr), the inclination distributions obtained in the model 
end up being incorrect. See, for example, the case with aN,o = 28 AU and r = 30 Myr in 
Figure inb. This specific case indicates cxj = 6-10°; it can be ruled out at a 99.8% conhdence 
level from the CFEPS detections. 

For aN,o = 24 AU, on the other hand, the 2:1 and 7:4 resonances are at 38.1 and 34.9 AU. 
These resonances therefore cannot deposit bodies into the main belt during the initial stages 
of migration (because they are not located in the main belt during these initial stages). This 
means that the orbital inclinations in the ISR region can be excited by Neptune’s encounters 
before the main phase of the implantation starts. It starts when the 2:1 resonance moves 
beyond ~40 AU, or equivalently, when Neptune moves beyond a ~ 40/2^/^ = 25.2 AU. 
When exactly this happens depends both on on.o and r. For example, with aN,o = 24 AU 
and r = 30 Myr, Neptune moves past 25.2 AU at f ~ 10 Myr after the start of migration. 
This delay is long enough to raise the mean orbital inclination in the ISR region to ~10°. 
This explains why the orbital inclinations of HCs are generally higher when aN,o is smaller. 

Figure [10] shows the eccentricity (panel c) and inclination (panel d) distributions of 
HCs obtained for aN,o = 24 AU and r = 10 Myr. While the eccentricity distribution is 
(slightly) discrepant when compared to the CFEPS detections, the inclination distribution 
looks good (K-S probability 29%). Unlike in Figure |5ji, panel d of Figure [TOl compares the 
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inclination distribntions all the way down to 5°. The CFEPS inclination distribntion is steep 
from 5° to 10°, and shallow above 10°. This may suggest that the underlying distribution 
is bimodal, perhaps because it resulted from capture at two different stages of Neptune’s 
migration (see discussion in Section 6). The model distribution obtained with aN,o = 24 AU 
and r = 10 Myr is a good proxy for the overall shape of the CFEPS curve. We conclude 
from the inclination distribution of HCs that Neptune’s migration was slow (r > 10 Myr) 
and long range (on.o ^ 25 AU). 

Figures [TTh shows the mean orbital inclination of bodies in the ISR as a function of 
time. The mean inclination of the ISR population steadily increases with time. It is ~5° at 
10® yr, ~10° at 10^ yr, ~15° at 3 x 10^ yr, and ~20° at 10® yr. This makes it obvious that 
bodies captured in the main belt late will have, on average, larger orbital inclinations than 
bodies captured early, and explains the trends discussed above (see also Figure [T^ . Note, 
however, that the number of bodies available in the ISR drops for f > 10 Myr fFigure fTTbb 
This is because bodies in the scattered disk evolve to very long orbital periods, or move to 
short orbital periods and are subsequently ejected from the solar system by Jupiter. This 
means that the number of bodies available for capture at very late times is relatively small. 
The very late captures {t > 10® yr) are therefore not very important for the overall statistics. 

In summary, we hnd that the inclination distribution of the HCs is a reflection of the 
inclination distribution in the ISR, weighted by the number of bodies in the ISR, and time 
integrated over the capture window that depends both on on.o and r. In addition, the 
inclination distribution becomes modihed by the dynamical processes involved in step 2 of 
the capture process (see Section 5.2). Our main conclusion from the inclination constraint is 
that Neptune’s migration was long-range (on.o ^ 25 AU), and that the migration timescale 
was long (r > 10 Myr). The case with aN,o = 24 AU and r = 10 Myr works relatively 
well for the HCs, but it fails for Plutinos (Figure [TU]) . The case with aN,o = 24 AU and 
r = 30 Myr works for HCs with i > 10 deg and Plutinos. As we discuss in Section 6, 
the somewhat different timescales indicated by the HCs and Plutinos may be related to a 
two-stage migration of Neptune, with faster migration (r ~ 10 Myr) during the hrst stage 
and slower migration (r ~ 30 Myr) during the second. 


5.4. The Implantation Efficiency 

The efficiency of implantation of the disk bodies into the Kuiper belt is a product 
of partial efficiencies of the three steps described in Section 5.2. During the first step, 
bodies from the disk at < 30 AU are scattered by Neptune to > 30 AU, where they can be 
captured into resonances. The number of scattered bodies available for capture is a function 
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of time fFigureITTb). The capture efficiency in a resonance mainly depends on the resonance 
strength (e.g., the strong 3:2 resonance is expected to capture more bodies) and Neptune’s 
migration speed. In step 2, the resonant bodies can evolve to orbits with lower eccentricities 
assuming there is enough time for the secular cycles to act. This depends on how Neptune’s 
migration timescale compares with the period of secular cycles in a specific resonance. Also, 
since the secular cycles depend on Neptune’s eccentricity, this stage is affected by Neptune’s 
eccentricity behavior during migration. Finally, whether a resonant orbit is or is not released 
from a resonance during step 3, is mainly influenced by Neptune’s migration speed (more 
bodies are released for higher speeds). Moreover, the implantation of bodies into the main 
belt occurs only when the relevant resonances are present in the 40 < a < 47 AU region. 
For implantation via the 2:1 resonance, this requires that Neptune is beyond ~25 AU. 

Some of the trends can be identified in our results. For example, with aN,o = 24 AU, 
the efficiency of implantation on a stable orbit in the 3:2 resonance is P 3-,2 = 9.2 x 10 ^ for 
r = 10 Myr, = 5.3 x 10““^ for r = 30 Myr, and P 3:2 = 2.0 x 10““^ for r = 100 Myr. Also, 
for r = 30 Myr, P^.^ = 5.3 x lO”"^ for on.o = 24 AU, P 3-,2 = 1.2 x 10“^ for aN,o = 26 AU, and 
F 3:2 = 2.3 X 10“^ for aN,o = 28 AU. Thus, a longer migration timescale leads to lower P 3 . 2 , 
and larger aN,o leads to higher P 3 . 2 . The trends for the implantation in the HC region are 
similar. For example, with r = 30 Myr, Phc = 1-9 x lO”"^ for on.o = 24 AU, Pro = 2.1 x lO”'^ 
for aN,o = 26 AU, and Prc = 1.1 x 10 “^ for on.o = 28 AU@ 

In our preferred case (aN,o = 24 AU and r = 10 or 30 Myr), the main-belt capture 
efficiency is Pro — 2-4 x 10“^ for each initial particle in the original disk. With M^isk = 
20 MEarth, the total mass of the hot population would therefore be Mrc = 0.004-0.008 MEarth- 
This is satisfactory when compared to Mrc — 0.01 Me estimated by Fraser et ah (2014), 
especially because Fraser’s estimate has a considerable uncertainty. Also, scaling from Jupiter 
Trojans, there should have been (3-4) x 10^ planetesimals with absolute magnitude H < 9 
in the original disk (Nesvorny et al. 2013). With Prc — 2-4 x 10“^, our model would predict 
~7000-14,000 HCs with P < 9, while Adams et al. (2014) give 19, 000 ± 5, 000 from the 
Deep Ecliptic Survey (DES) for the whole population of the classical belt. This agreement 
is reasonable. [We note that the model estimates discussed here were obtained with cn.o = 0 
and whenever cn stayed low during the migration. The cases with eN,o = 0.1 and/or the 
ones where the mean motion resonances between Uranus and Neptune acted to temporarly 
increase cn during the migration tend to produce larger Prc, by a factor of a few, but more 
simulations would need to be done to establish this trend convincingly. The trend could be 


®p 3:2 (or Pro) is the probability that an original disk object ends up on a stable orbit in the 3:2 resonance 
(or in the HC region). We compute this probability by dividing the number of 3:2 resonant (or HC) bodies 
at the end of our simulations (4 Gyr) by the number of bodies in the original disk (10®). 
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related to the dependence of the secnlar cycles inside the mean motion resonances on cn.o-] 

A major problem is identihed when we consider the captnre probability of the resonant 
objects. With P 3.2 ~ 5 x 10“^ for the preferred case with aN,o = 24 AU and r = 30 Myr, we 
would predict that the 3:2 resonance population should host ~ 1.5-3 times more objects than 
the hot population in the main belt, while according to the CFEPS survey there are ~3.5 
times as many HCs as Plutinos (with absolute magnitude H < 8 ; B. Gladman, personal 
communication). This would indicate P 3 : 2 /Phc — 0.3 (this estimate has a <50% formal 
uncertainty). The 3:2 resonance is thus obviously overpopulated in our simulations, roughly 
by a factor of 5-10. We call this the resonance overpopulation problem. This problem was 
already noted in many previous dynamical models of Kuiper belt formation (e.g., Hahn & 
Malhotra 2005, Levison et ah 2008, Morbidelli et al. 2008). 

There are several potential solutions to this problem. For example, we performed several 
simulations with r = 100 Myr and found that Pro tends to be higher, and P 3:2 tends to be 
lower, than in the cases with r = 30 My. For example, with aN,o = 26 AU and r = 100 Myr 
(see Figure US] for the orbital distribution of bodies obtained in this simulation), we found 
that Pro — 9 x 10 “"^ and = 1.6 x 10 “^, thus indicating P^-, 2 /Pnc — 5.6. Also, with 
aN,o = 28 AU and r = 100 Myr, Prc ^ 3.2 x 10“^ and P^a = 1-4 x 10“^, so Pz-. 2 /Pnc — 2.3. 
These ratios are more similar to the value P 3 : 2 /Phc ~ 3.5 inferred from observations. The 
case with aN,o = 26 AU produces a slightly wider inclination distribution than indicated 
by observations, while the one with aN,o = 26 AU produces a much narrower inclination 
distribution, thus suggesting the possibility that an intermediate value of on.o — 26-27 AU 
would give the correct result. We do not give much emphasis to the cases with r = 100 Myr in 
this paper, because it is not clear how these cases relate to our current models of the planetary 
instability and migration (e.g., they may require a very low mass of the planetesimal disk). 
New planetary instability simulations, with a focus on very slow migration timescales, will 
need to be performed. Another solution of the resonance overpopulation problem, which will 
also require additional modeling effort that goes beyond the scope of this paper, is discussed 
in Section 6 . 

We point out that the overpopulation problem is not specihc to the 3:2 resonance. In¬ 
stead, nearly all resonances are overpopulated. On the other hand, when we compare, relative 
to each other, the number of bodies captured in different resonances in our simulations, we 
hnd that these populations have roughly the right proportions. For example, for on.o = 24 
AU and r = 30 Myr, we find that Pz: 2 /P 2 :i — 3 and Pz: 2 /P 2 :i — 7. This is comparable to 
the resonance population statistics in the 3:2, 2:1 and 5:2 resonances discussed in Gladman 
et al. (2014) (even though the population in the 5:2 resonance was previously thought to be 
comparable to that in the 3:2 resonance; Gladman et al. 2012). 
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Another notable result obtained from our simulations concerns the NTs. Two problems 
were identihed in previous modeling of the NT capture: (1) the capture efficiency obtained 
in the previous simulations was orders of magnitude too high, and (2) the inclination 
distribution was too narrow (e.g., Nesvorny & Vokrouhlicky 2009, Parker 2015). Related to 
(1), only hve NTs were captured out of the original 10® disk particles for aN,o = 24 AU and 
r = 30 Myr. This indicates the capture probability Pnt ~ 5 x 10“® and is more in line with 
observations (Alexandersen et ah 2014). As for (2), four of hve stable NTs produced by the 
reference run have inclinations > 20° (Figure 0]). This is encouraging, but better statistics 
will be needed to compare things more carefully. 


6. Discussion and Conclusions 

The Gomes mechanism was identihed here to have fundamental importance for the origin 
of dynamical structure in the Kuiper belt. The basic requirement for the Gomes mechanism 
to work is that the migration timescale r is comparable to, or longer than, the secular cycles 
inside the mean motion resonances such as 2:1, 3:2, 7:3, etc. With on.o 25 AU, the 2:1 
resonance is initially below 40 AU, sweeps over the main belt location at 40 < a < 47 AU 
during Neptune’s migration, and is responsible for the delivery of most objects into the main 
belt. Since the secular cycles in the 2:1 resonance have a several-Myr period, r needs to be at 
least several Myr for this to work. If, instead, r ~ 1 Myr, objects can still be captured in the 
main belt region, assuming that Neptune had an orbit with substantial orbital eccentricity 
(Levison et ah 2008). In this case, orbits evolve from the scattered disk to the main belt by 
normal secular cycles outside the mean motion resonances (Dawson & Murray-Glay 2012). 
These cycles have a shorter period and work for shorter migration timescales. 

A major problem with the high-eccentricity (cn > 0.1) phase of Neptune is the opposing 
constraints from the hot and cold populations, as explained in Dawson & Murray-Glay 
(2012). Specihcally, to preserve the GGs, the eccentricity of Neptune cannot be large and/or 
must be damped fast (Batygin et al. 2011, Wolff et ah 2012). To capture the HGs by 
the normal secular cycles outside the mean motion resonances, however, Neptune’s initial 
eccentricity must be relatively high at some point during Neptune’s migration. While these 
two constraints rule out most of parameter space, Dawson & Murray-Glay (2012), working 
under the assumption that the HGs were captured by the normal secular cycles, found 
solutions that satisfy both. This niche of parameter space is eN,o — 0.1 and on.o > 28 AU, 
essentially meaning that Neptune’s migration would have to be short range. 

The beauty of the Gomes mechanism is that it works even if Neptune’s eccentricity was 
never large. This is because this mechanism does not rely on the normal secular eccentric- 
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ity oscillations forced by eccentric Neptnne on orbits ontside the mean motion resonances. 
Instead, it appears as a prodnct of large eccentricity oscillations dne to the existence of 
large-amplitnde secnlar cycles inside the mean motion resonances, akin to those hrst pointed 
ont for the 3:1 Jnpiter resonance in the asteroid main belt (Wisdom 1982). Here we showed 
that the Gomes mechanism is the dominant implantation mechanism from < 30 AU in a 
regime where cn is low. Therefore, the argnment of Dawson & Mnrray-Clay (2012) that 
eN,o > 0.12 is reqnired to explain the HCs does not apply. On a related note, Levison et al. 
(2008) snggested that the existence of the detached disk with orbital perihelia extending to 
~40 AU is a conseqnence of an eccentric phase of Neptnne, for Neptnne to be capable of 
scattering objects to g ~ 40 AU, and nsed cn = 0.3 in their simnlations. Here we showed 
that the Gomes mechanism can prodnce the correct orbital architectnre of the detached disk 
even if cn stays low, assnming that the migration timescale was long (e.g., Fignres 4 and 
13). 

In this work, we stressed the importance of the inclination distribntion of the KBOs. 
This is becanse the inclination distribntion has been relatively well characterized from ob¬ 
servations and can therefore be nsed to constrain models. We showed that the Gomes mech¬ 
anism is capable of prodncing the observed wide inclination distribntion from a dynamically 
cold disk at <30 AU, assnming that Neptnne’s migration was long-range (aN,o ^ 25) AU and 
slow (r > 10 Myr)@ Since the Gomes mechanism is insensitive to Neptnne’s eccentricity, the 
eccentricity conld have been negligible dnring Neptnne’s migration in mnch the same way 
as originally proposed by Malhotra (1993, 1995) and later nsed by Hahn & Malhotra (1999, 
2005) to model the origin of the Kniper belt from a disk at > 30 AU. 

Does this mean that the Nice-type instability never happened? Not really. On one 
hand, the resnlts presented here seem to rnle ont the strong instability version of the Nice 
model, where Neptnne was thrown onto an eccentric orbit with a > 25 AU. This is becanse, 
as we discnssed above, the short-range migration of Neptnne into a dynamically cold disk 
at < 30 AU wonld lead to a narrow distribntion of orbital inclinations, in contradiction 
to observations (also see Levison et al. 2008). Moreover, to stabilize eccentric Neptnne at 
> 25 AU, the disk wonld have to be massive (~ 50 MEarth; Nesvorny & Morbidelli 2012; 
hereafter NM12), and wonld prodnce fast migration of Neptnne, thns leading to a donble 
contradiction, becanse both r and aN,o wonld be ont of the plansible range identihed here 
(Neptnne’s migration needs to be slow, not fast, to explain the inclination distribntion). 


®Using an initially strongly excited disk is counterproductive, because this has the consequence, as shown 
by our additional simulations of Neptune migrating into a pre-heated disk, that the implantation efficiency 
in the main belt drops by a factor of several. 
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On the other hand, some dynamical instability in the onter solar system clearly mnst 
have happened. The best evidence for this is the eccentric orbit of Jnpiter, which can 
be conveniently explained if Jnpiter snffered enconnters with an ice giant (Morbidelli et ah 
2009a). A discontinnons evolntion of Jnpiter’s semimajor axis, known as the jumping-Jupiter 
model, presnmably prodnced by varions scattering events dnring the epoch of planetary 
enconnters, is also reqnired from the terrestrial planet (Brasser et ah 2009, 2013; Agnor & 
Lin 2012) and asteroid belt constraints (e.g., Morbidelli et ah 2010). Moreover, planetary 
enconnters may be needed to explain the captnre and orbital distribntion of Jnpiter Trojans 
and irregnlar satellites (Nesvorny et ah 2007, 2013, 2014). 

A new model of planetary instability has recently been proposed (Nesvorny 2011, NM12, 
Batygin et al. 2012). This model has been the framework of several newer pnblications 
discnssed above. It postnlates that the early solar system had an extra ice giant, which was 
ejected into interstellar space dnring the instability. Fignre [H] illnstrates this possibility. 
In this model, five onter planets start in a relatively relaxed confignration with Neptnne at 
~ 22 AU. The first thing that happens in the simnlation is that Neptnne migrates into the 
onter disk located at 24-30 AU. After Neptnne reaches ~ 28 AU, the instability happens, 
dnring which the extra ice giant has enconnters with all the other onter planets, and is 
snbseqnently ejected by Jnpiter. The main featnres of this model which are most relevant 
for the Kniper belt are: (1) Neptnne’s eccentricity and inclination are never large (e < 0.1 
and i < 2°), (2) the initial mass of the onter disk at < 30 AU is relatively small (~ 15- 
20 MEarth; NM12), therefore implying a slow migration of Neptnne, (3) Neptnne’s semimajor 
axis discontinnonsly changes (by ~0.2-0.5 AU) when Neptnne is at ~28 AU, as a resnlt of 
one or two very close enconnters with the ejected ice giant; the migration rate is slower 
after the jnmp than it was dnring the previons migration stage, and (4) the ejected ice giant 
briefly overlaps with the Kniper belt (Batygin et al. 2012). 

As for (4), we carefnlly looked into several instability cases from NM12 and fonnd that 
the relevant period dnring which the ice giant’s orbit overlaps with the Kniper belt is too brief 
to signihcantly affect the orbits in the Kniper belt. Batygin et al. (2012), who fonnd larger 
effects in abont 50% of stndied cases, did so probably becanse their work covered a broad 
range of possibilities, with some of their instability cases being somewhat too cataclysmic, 
in onr opinion, to prodnce the solar system as we know it now. Items (1) and (2) present 
the right conditions for the Gomes mechanism to work and play a dominant role over other 
implantation mechanisms. 

Unlike in the idealized case stndied here, Neptnne’s migration in Fignre [HI and other 
cases reported in NM12, happens in two stages. Dnring the hrst stage, that is, before the 
instability happens, Neptnne migrates with r ~ 10 Myr for Mdisk = 20 MEarth or r ~ 20 Myr 
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for Mdisk = 15 MEarth (these are the best exponential fits in the NM12 cases we looked at). 
During the second stage, that is after the instability, Neptune migrates with r ~ 30 Myr for 
Md isk = 20 Msarth or r ~ 50 Myr for Md^k = 15 MEarth- These best-£t r values are only 
approximate, because the real migration is not exactly exponential, and the effective r is 
typically longer as time progresses. 

These timescales, and the long-range nature of Neptune’s migration in NM12, agree 
quite nicely with the constraints on Neptune’s migration derived from the Kuiper belt in 
this work. This shows that the NM12 instability model, which was developed entirely from 
constraints unrelated to the Kuiper belt, may have some relevance for the early evolution of 
the solar system. 

The inclination distribution of the HCs may provide some evidence for the two stage 
migration of Neptune in the NM12 model. This is because the HCs can be captured into the 
main belt both during the hrst phase, when Neptune’s migration was faster, and during the 
second phase, when the migration was slower. The main implication of this is that the HCs 
can be a composite of two populations captured at two different stages. From the discus¬ 
sion in this paper, these populations are expected to have different inclination distributions, 
thus potentially explaining why the CFEPS detections of the HCs show different slopes for 
5° < i < 10° (captures during the hrst stage) and i > 10° (second stage). Also, Plutinos and 
other resonant populations, captured during the hrst stage, would be released when Nep¬ 
tune jumped during the NM12 instability. This could relieve the resonance overpopulation 
problem discussed in Section 5.4. The present resonant populations would then have to be 
captured entirely during the second phase, when the migration of Neptune was slower. This 
could explain why we are seeing a preference in our results for slightly longer r values for 
Plutinos than for HCs (e.g., the case with r = 10 Myr is works well for the HCs but does 
not really work for Plutinos; Figure [TUD . 

Much work has yet to be done to fully understand the dynamics of the KBOs during 
Neptune’s migration, and how the history of Neptune’s orbit is constrained by the dynamical 
structure of the Kuiper belt. Clearly, the idealized migration model studied here is a major 
simplification. We used it to highlight several interesting results that can be obtained within 
the framework of this model. We now plan to increase the realism of the model by considering 
the two-stage migration from NM12. We believe that this will be a crucial step toward 
resolving the resonance overpopulation problem that plagued previous studies of Kuiper 
belt formation. In a companion paper, Nesvorny (2015), we study the implications of the 
NM12 instability model for the cold classical belt. 


This work was supported by NASA’s Outer Planet Research (OPR) program. All CPU- 


intensive simulations in this work were performed on NASA’s Pleiades Supercomputer!^ We 
thank W. F. Bottke, L. Dones, B. Gladman, H. F. Levison, A. Morbidelli, D. Vokrouhlicky, 
and an anonymous reviewer for helpful comments on this work. 


REFERENCES 

Adams, E. R., Gulbis, A. A. S., Elliot, J. L., et ah 2014, AJ, 148, 55 
Agnor, G. B., & Lin, D. N. G. 2012, ApJ, 745, 143 

Alexandersen, M., Gladman, B., Kavelaars, J. J., et ah 2014, arXiv:1411.7953 
Allen, R. L., Gladman, B., Kavelaars, J. J., et ah 2006, ApJ, 640, L83 
Batygin, K., Brown, M. E., & Fraser, W. G. 2011, ApJ, 738, 13 
Batygin, K., Brown, M. E., & Betts, H. 2012, ApJ, 744, L3 
Bernstein, G. M., Trilling, D. E., Allen, R. L., et ah 2004, AJ, 128, 1364 
Brasil, P. I. O., Gomes, R. S., & Soares, J. S. 2014a, A&A, 564, 44 
Brasil, P. I. O., Nesvorny, D., & Gomes, R. S. 2014b, AJ, 148, 56 

Brasser, R., Morbidelli, A., Gomes, R., Tsiganis, K., & Levison, H. F. 2009, A&A, 507, 1053 
Brasser, R., Walsh, K. J., & Nesvorny, D. 2013, MNRAS, 433, 3417 
Brown, M. E. 2001, AJ, 121, 2804 

Brown, M. E., Schaller, E. L., & Fraser, W. G. 2011, ApJ, 739, 60 

Brucker, M. J., Grundy, W. M., Stansberry, J. A., et ah 2009, Icarus, 201, 284 

Ghiang, E. h, & Jordan, A. B. 2002, AJ, 124, 3430 

Ghiang, E. h, Jordan, A. B., Millis, R. L., et ah 2003, AJ, 126, 430 

Dawson, R. h, & Murray-Glay, R. 2012, ApJ, 750, 43 

Duncan, M. J., Levison, H. F., & Budd, S. M. 1995, AJ, 110, 3073 


^°http://www.nas.nasa.gov/hecc/resources/pleiades.htinl 



Fraser, W. C., Brown, M. E., Morbidelli, A., Parker, A., & Batygin, K. 2014, ApJ, 782, 100 

Gladman, B., Marsden, B. G., & Vanlaerhoven, G. 2008, The Solar System Beyond Neptune, 
43 

Gladman, B., Lawler, S. M., Petit, J.-M., et al. 2012, AJ, 144, 23 

Gladman, B., Bannister, M., Kavelaars, J., et al. 2014, AAS/Division for Planetary Sciences 
Meeting Abstracts, 46, #507.01 

Gomes, R. S. 2003, Icarus, 161, 404 

Gomes, R. S. 2011, Icarus, 215, 661 

Gomes, R. S., Morbidelli, A., & Levison, H. F. 2004, Icarus, 170, 492 

Gomes, R. S., Gallardo, T., Fernandez, J. A., & Brunini, A. 2005, Gelestial Mechanics and 
Dynamical Astronomy, 91, 109 

Gulbis, A. A. S., Elliot, J. L., Adams, E. R., et al. 2010, AJ, 140, 350 
Hahn, J. M., & Malhotra, R. 1999, AJ, 117, 3041 
Hahn, J. M., & Malhotra, R. 2005, AJ, 130, 2392 

Kavelaars, J., Jones, L., Gladman, B., Parker, J. W., & Petit, J.-M. 2008, The Solar System 
Beyond Neptune, 59 

Kavelaars, J. J., Jones, R. L., Gladman, B. J., et al. 2009, AJ, 137, 4917 

Kenyon, S. J., Bromley, B. G., O’Brien, D. P., & Davis, D. R. 2008, The Solar System 
Beyond Neptune, 293 

Knezevic, Z., Milani, A., Farinella, P., Froeschle, G., & Froeschle, G. 1991, Icarus, 93, 316 
Kozai, Y. 1962, AJ, 67, 591 

Kuchner, M. J., Brown, M. E., & Holman, M. 2002, AJ, 124, 1221 
Levison, H. F., & Duncan, M. J. 1994, Icarus, 108, 18 
Levison, H. F., & Stern, S. A. 2001, AJ, 121, 1730 
Levison, H. F., & Morbidelli, A. 2003, Nature, 426, 419 



- 26 - 

Levison, H. F., Morbidelli, A., Vanlaerhoven, C., Gomes, R., & Tsiganis, K. 2008, Icarus, 
196, 258 

Li, J., Zhou, L.-Y., & Sun, Y.-S. 2014, MNRAS, 437, 215 
Lykawka, P. S., & Mukai, T. 2008, AJ, 135, 1161 
Malhotra, R. 1993, Nature, 365, 819 
Malhotra, R. 1995, AJ, 110, 420 

Morbidelli, A., Thomas, F., & Moons, M. 1995, Icarus, 118, 322 

Morbidelli, A., Tsiganis, K., Crida, A., Levison, H. F., & Gomes, R. 2007, AJ, 134, 1790 

Morbidelli, A., Levison, H. F., & Gomes, R. 2008, The Solar System Beyond Neptune, 275 

Morbidelli, A., Brasser, R., Tsiganis, K., Gomes, R., & Levison, H. F. 2009a, A&A, 507, 
1041 

Morbidelli, A., Levison, H. F., Bottke, W. F., Dones, L., & Nesvorny, D. 2009b, Icarus, 202, 
310 

Morbidelli, A., Brasser, R., Gomes, R., Levison, H. F., & Tsiganis, K. 2010, AJ, 140, 1391 

Morbidelli, A., Gaspar, H. S., & Nesvorny, D. 2014, Icarus, 232, 81 

Murray-Glay, R. A., & Ghiang, E. 1. 2005, ApJ, 619, 623 

Murray-Glay, R. A., & Ghiang, E. 1. 2006, ApJ, 651, 1194 

Nesvorny, D. 2011, ApJ, 742, L22 

Nesvorny, D., & Morbidelli, A. 2012 (NM12), AJ, 144, 117 

Nesvorny, D., & Roig, F. 2000, Icarus, 148, 282 

Nesvorny, D., & Roig, F. 2001, Icarus, 150, 104 

Nesvorny, D., Vokrouhlicky, D., & Morbidelli, A. 2007, AJ, 133, 1962 

Nesvorny, D., Vokrouhlicky, D., & Morbidelli, A. 2013, ApJ, 768, 45 

Nesvorny, D., Vokrouhlicky, D., & Deienno, R. 2014, ApJ, 784, 22 

Nesvorny, D. 2015. Submitted to AJ 



-27- 

Noll, K. S., Grundy, W. M., Chiang, E. L, Margot, J.-L., & Kern, S. D. 2008a, The Solar 
System Beyond Neptune, 345 

Noll, K. S., Grundy, W. M., Stephens, D. G., Levison, H. F., & Kern, S. D. 2008b, Icarus, 
194, 758 

Parker, A. H. 2015, Icarus, 247, 112 

Parker, A. H., & Kavelaars, J. J. 2010, ApJ, 722, L204 

Petit, J.-M., Kavelaars, J. J., Gladman, B. J., et ah 2011, AJ, 142, 131 

Petit, J.-M., Kavelaars, J. J., Gladman, B., Jones, L., & Parker, J. 2014, AAS/Division for 
Planetary Sciences Meeting Abstracts, 46, ^^507.07 

Petit, J.-M. et al. 2015. To be snbmitted. 

Press, W. H., Tenkolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Gambridge: 
University Press, 2nd ed. 

Stern, S. A., & Golwell, J. E. 1997, AJ, 114, 841 

Tegler, S. G., & Romanishin, W. 2000, Nature, 407, 979 

Thommes, E. W., Duncan, M. J., & Levison, H. F. 1999, Nature, 402, 635 

Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. F. 2005, Nature, 435, 459 

Wisdom, J. 1982, AJ, 87, 577 

Volk, K., & Malhotra, R. 2011, ApJ, 736, 11 

Wolff, S., Dawson, R. L, & Murray-Glay, R. A. 2012, ApJ, 746, 171 


This preprint was prepared with the AAS IATJtiX macros v5.0. 




30 40 50 60 

Semimajor axis (AU) 

Fig. 1.— The orbital elements of KBOs observed in three or more oppositions. Various 
dynamical classes are highlighted. The HCs with i > 5° and NTs are denoted by larger 
dots, and the CCs are denoted by smaller dots. Note the wide inclination distribution of the 
HCs in panel (b) with inclinations reaching above ~ 30°. The solid lines in panel (a) follow 
the borders of important mean motion resonances. For the NTs, we show the approximate 
location of stable librations from Nesvorny & Dones (2002). The low-inclination orbits with 
40 < a < 42 AU are unstable due to secular resonance overlap (z/y and z/g; Knezevic 1991, 
Duncan et al. 1995). 






Inclination (deg) 

Fig. 2.— The cumulative inclination distribution of Plutinos in the 3:2 resonance with 
Neptune. The distribution of 29 Plutinos detected by CFEPS (dashed line, Petit et ah 2011) 
is compared to a model distribution (solid line). The K-S test applied to these distributions 
shows that the likelihood that they can be obtained from the same underlying distribution 
is 2 X 10“^^. This rules out the model. The model distribution is a result of a numerical 
simulation where we considered Neptune’s migration into a dynamically cold disk (ue = 0.1 
and (Tj = 2°) at a < 30 AU. See Section 3 for a description of the model. Here we used 
aN,o = 28 AU, eN,o = 0.1 and iN,o = 0.67°. The radial migration and eccentricity damping 
were applied to Neptune’s orbit on an e-folding timescale r = 1 Myr. To compare apples with 
apples, the CFEPS simulator (see Section 4) was used to compute the detection statistics 
from the population of bodies that survived in the 3:2 resonance at the end of the simulation. 
Thus, both distributions shown in this plot include the observational bias of CFEPS. The 
inclination distribution obtained for larger values of eN,o (Levison et al. 2008 used eN,o = 0.3) 
is similar to the one shown here for eN,o = 0.1. 
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Fig. 3.— The same as Figured] but without labeling of different populations. This plot is 
useful for a visual comparison with the model results. 
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Fig. 4.— The orbital elements of bodies captured in the Kuiper belt in a model with 
aN,o = 24 AU and r = 30 Myr. The HCs and NTs are denoted by larger symbols. 
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Fig. 5.— The cumulative distribution of eccentricities (left) and inclinations (right) for 
Plutinos (upper) and HCs (lower). The dashed lines show the actual CFEPS detections (29 
Plutinos, and 10 HCs with i > 10°). The solid lines show the distributions of model bodies 
(on,o = 24 AU and r = 30 Myr) detected by the CFEPS simulator. Both the observed and 
model distributions plotted here therefore contain the CFEPS observational bias. For the 
HCs, we compare the distributions for i > 10° to avoid any potential contamination of the 
detection statistics from the CCs, which are not modeled here. 
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Fig. 6.— An illustration of the Gomes implantation mechanism. The panels show: the (a) 
path of a disk particle in the (a, e) projection; the two red dots show the initial and hnal 
orbits, (b) semimajor axis, (c) eccentricity, (d) inclination, (e) perihelion argument cn, and (f) 
2:1 resonant angle a 2 :i = 2A —An — ci7, where A and An are the particle’s and Neptune’s mean 
longitudes, and vu is the particle’s perihelion longitude. After being scattered by Neptune 
and experiencing Kozai cycles inside the 2:1 mean motion resonance, the disk particle ends 
up on a high-inclination orbit in the main belt. 
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Fig. 7.— Capture of a disk particle on a stable, high-inclination orbit in the 3:2 resonance. 
The panels show: the (a) path of the disk particle in the (a, e) projection; the two red dots 
show the initial and hnal orbits, (b) semimajor axis, (c) eccentricity, (d) inclination, (e) 
perihelion argument u, and (f) 3:2 resonant angle a 3-,2 = 3A — 2An — cu. 
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Fig. 8.— Capture of a disk particle on a high-inclination Neptune Trojan orbit. The panels 
show: the (a) path of the disk particle in the (a, e) projection; the two red dots show the 
initial and final orbits, (b) semimajor axis, (c) eccentricity, (d) inclination, (e) perihelion 
argument u, and (f) 1:1 resonant angle ai,! = A — An- 
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Fig. 9.— The dependence of the inclination distribution on r and aN,o: (a) Plutinos, and (b) 
HCs. These distributions are the intrinsic distributions obtained in the model (i.e., they do 
not include any observational bias). The inclination distribution of Plutinos in panel (a) is 
sensitive to the assumed migration timescale. The labels in panel (a) denote the distributions 
obtained for different values of r. In panel (b), we show the inclination distribution of HCs 
obtained with r = 30 Myr. The two lines correspond to different starting positions of 
Neptune (aN,o = 24-26 AU and 28 AU; the cases with aN,o = 24 and 26 AU were similar and 
were put together in this plot). 
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Fig. 10.— The cumulative distributions of eccentricities (left) and inclinations (right) for 
Plntinos (npper) and HCs (lower). The dashed lines show the actnal CFEPS detections (29 
Plntinos, and 21 HCs with i > 5°). The solid lines show the distribntions of model bodies 
(on,o = 24 AU and r = 10 Myr) detected by the CFEPS simulator. Both the observed 
and model distribntions plotted here therefore contain the observational bias of CFEPS. For 
the HCs we compare the distribntions for i > 5° to avoid contamination of the detection 
statistics from the CCs, which are not modeled here. 
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Fig. 11.— The mean orbital inclination and number of objects in the Intermediate Source 
Region (ISR), dehned as 40 < a < 47 AU and q < Qn- The lower panel shows that it 
takes 10® to 10^ yr to build up the ISR population, which then decays by about an order of 
magnitude over the next 10® years. The mean orbital inclination of the ISR population in 
panel (b) steadily increases with time. The different lines in panels (a) and (b) correspond 
to models with different r and aN,o- The overall shape of the lines is insensitive to these 
parameters. 
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Fig. 12.— The inclination distribution of objects in the IRS for t = 1, 10, 30, 100, and 300 
Myr after the start of Neptune’s migration. The plot shows how the orbital inclinations in 
the ISR increase as a result of encounters with Neptune. Here we used aN,o = 24 AU and 
r = 30 Myr. 
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Fig. 13.— The orbital elements of bodies captured in the Kuiper belt in a model with 
flN,o = 26 AU and r = 100 Myr. The HCs and NTs are denoted by larger symbols. Note 
that the number of bodies captured in the HC region is ~5 times larger than in Figure 01 The 
detached disk beyond 50 AU is also more populated and extends to g > 40 AU. Both these 
results are a consequence of the very slow migration of Neptune assumed in this simulation. 
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Fig. 14.— The orbital histories of the giant planets from the instability/migration model of 
NM12. In this example, the hfth giant planet was initially placed on an orbit between Satnrn 
and Uranus and was given a mass equal to Neptune’s mass. The three inner planets were 
started in the (3:2,4:3) resonant chain, and Uranus and Neptune on non-resonant orbits with 
ou = 17 AU and on = 22 AU. Ten thousand particles, representing the outer planetesimal 
disk, were distributed with semimajor axes 23.5 < a < 29 AU, surface density S = 1/a, 
and low eccentricities and low inclinations. With the total disk mass Maisk = 20 Me, each 
disk particle has approximately one Pluto mass, (a) The semimajor axes (solid lines), and 
perihelion and aphelion distances (thin dashed lines) of each planet’s orbit. The inner ice 
giant was ejected into interstellar space at t ~ 10 Myr after the start of the simulation. 
The hnal orbits of the four remaining planets are a good match to those in the present 
solar system (thin dashed lines), (b) The period ratio Pn/Pu- The thin dashed line shows 
Tn/Pu = 1-96 corresponding to the present orbits of Uranus and Neptune. Panels (c) and 
(d) show the eccentricity and inclination of Neptune’s orbit. The bold dashed lines in panels 
(a) and (c) approximate Neptune’s migration with Tq = 30 Myr and damping of Neptune’s 
eccentricity with Tg = 7 Myr. These long migration/damping timescales are characteristic 
for the migration/instability models from NM12. 




























